interaction_plot <- function(model,var1,var2,int,vcov,ci=.95,
                   xlab=var2,ylab=paste("Marginal Effect of",var1),
                   main="",
                   me_lty=1,me_lwd=.75,me_col="black",
                   ci_lty=1,ci_lwd=.5,ci_col="black",
                   yint_lty=2,yint_lwd=.5,yint_col="black"){
  require(ggplot2)
  alpha <- 1- .95
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov(model)
  z0 <- seq(min(model.frame(model)[,var2],na.rm=T),max(model.frame(model)[,var2],na.rm=T),length.out=100)
  dy.dx <- beta.hat[var1] + beta.hat[int]*z0
  se.dy.dx <- sqrt(cov[var1,var1] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov[var1,ncol(cov)])
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  
  ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
    labs(x=xlab,y=ylab,title=main) +
    geom_line(aes(z0, dy.dx+2),size = me_lwd, 
              linetype = me_lty, 
              color = me_col) +
    geom_ribbon(aes(x = z0, ymin = lwr+2, ymax = upr+2), fill = "grey70", alpha = 0.5) +
    geom_hline(yintercept=2,linetype=yint_lty,
               size=yint_lwd,
               color=yint_col)
  
}




interaction_plot_binary <- function(model,var1,var2,int,vcov,ci=.95,
                             xlab=var2,ylab=paste("Marginal Effect of",var1),
                             main="",
                             me_lty=1,me_lwd=2,me_col="black",
                             ci_lty=1,ci_lwd=.5,ci_col="black",
                             yint_lty=2,yint_lwd=.5,yint_col="black"){
  require(ggplot2)
  alpha <- 1- .95
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov(model)
  z0 <- seq(min(model.frame(model)[,var2],na.rm=T),max(model.frame(model)[,var2],na.rm=T),length.out=2)
  dy.dx <- beta.hat[var1] + beta.hat[int]*z0
  se.dy.dx <- sqrt(cov[var1,var1] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov[var1,ncol(cov)])
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  
  ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
    labs(x=xlab,y=ylab,title=main) +
    geom_point(aes(z0, dy.dx),size = me_lwd, 
              color = me_col) +
    geom_errorbar(aes(ymin = lwr, ymax =  upr), colour = "black", width = 0.1) + 
    geom_hline(yintercept=0,linetype=yint_lty,
               size=yint_lwd,
               color=yint_col)
  
}






